Objetivos

TODO
TODO
TODO

Datos

En este projecto usaremos el fichero de datos “scrubbed.csv” obtenido en https://www.kaggle.com/datasets/NUFORC/ufo-sightings. Elimininaremos datos de tipo NA y cambiaremos paises vacios (i.e. pais == ““) por”Otros”. Tambien arreglamos el formato de los datos con la libreria lubridate. Usamos la libreria dplyr para eliminar columnas no necesarias.

# Leemos datos de CSV
datos <- read.csv("scrubbed.csv")

# Reemplacamos paises "" por "Otros"
datos$country <- ifelse(datos$country=="","otros",datos$country)

# Formateamos datos
datos$datetime <- as.POSIXct(datos$datetime, format =  "%m/%d/%Y %H:%M")
datos$'duration (seconds)' <- as.integer(datos$'duration..seconds.')
datos$latitude <- as.numeric(datos$latitude)

# Eliminamos datos NA
datos <- na.omit(datos)

library(dplyr)

datos <- datos %>% select(-one_of('duration..seconds.', 'duration..hours.min.', 'date.posted'))

Analisando los datos por pais, lógicamente podemos ver que la mayoria de los OVNIs reportados a la National UFO Reporting Center de EEUU ocurririon dentro de EEUU. Usaremos la libreria ggplot2 para generar un diagrama representando el numero de OVNIs reportados por pais y la libreria leaflet para generar un mapa con las longitudes y latitudes.

library(leaflet)
library(ggplot2)

by_country <- aggregate(datos$country, by=list(Country = datos$country), FUN=length)

ggplot(by_country, aes(x=by_country$Country, y=by_country$x, fill=by_country$Country)) +
  geom_bar(stat="identity", width=1) +
  xlab('Paises') +
  ylab('Avistamientos') +
  guides(fill = guide_legend(title = "Country" ))

leaflet(datos) %>%
  addTiles() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addMarkers(~longitude, ~latitude,
             popup = ~state, label = ~city,
             clusterOptions = markerClusterOptions())

Por eso, hemos decidido hacer nuestros analisis sobre datos solamente de EEUU. Tomamos una prueba de tamaño 1000.

datos_us <- datos[which(datos$country == 'us'),]

datos_us <- datos_us %>% select(-country)

datos_us <- datos_us[sample(nrow(datos_us), 1000), ]

Tema 3 - Estadística Descriptiva y Regresión

De la prueba de 1000 datos, se ha hecho el análisis solamente de las variables cuantitativas, las cuales son: (“datetime, latitude, longitude, duration(seconds)”).

Con el comando names(object…) se puede ver todas las columnas del fichero de datos y las referenciadas anteriormente.

names(datos_us)

## [1] "datetime"           "city"               "state"             
## [4] "shape"              "comments"           "latitude"          
## [7] "longitude"          "duration (seconds)"

Medidas de Posición - Tendencia Central y No Central

Con summary(object...) podemos ver la media, mediana y percentiles de varios datos, para calcular la moda usamos la funcion mode.

mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

datetime

summary(datos_us$datetime)

##                  Min.               1st Qu.                Median 
## "1947-04-15 14:00:00" "2001-11-03 17:36:15" "2007-07-24 05:37:30" 
##                  Mean               3rd Qu.                  Max. 
## "2004-10-25 21:59:05" "2011-09-09 05:12:30" "2014-05-06 22:30:00"

print(paste("Moda: ", mode(datos_us$datetime)))

## [1] "Moda:  2012-07-04 22:30:00"

latitude

summary(datos_us$latitude)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   19.73   34.05   38.74   38.22   41.85   61.22

print(paste("Moda: ", mode(datos_us$latitude)))

## [1] "Moda:  33.4483333"

longitude

summary(datos_us$longitude)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -158.19 -114.91  -91.64  -96.28  -81.03  -68.41

print(paste("Moda: ", mode(datos_us$longitude)))

## [1] "Moda:  -112.0733333"

duration (seconds)

summary(datos_us$'duration (seconds)')

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1      30     180    1086     600  109800

print(paste("Moda: ", mode(datos_us$'duration (seconds)')))

## [1] "Moda:  300"

Medidas de Dispersión

No se ha considerado necesario calcular la varianza y desviación poblacional

Medida	Formula
Cuasivarianza	var(x)
Cuasidesviación típica	sd(x)
Coeficiente de variación	coef_var(x)

coef_var <- function(x) {
  sd(x) / mean(x)
}

duration (seconds)

varianza <- var(datos_us$'duration (seconds)')
print(paste("Cuasivarianza: ", varianza))

## [1] "Cuasivarianza:  35067587.5226226"

sd <- sd(datos_us$'duration (seconds)')
print(paste("Cuasidesviación típica: ", sd))

## [1] "Cuasidesviación típica:  5921.78921632834"

cv <- coef_var(datos_us$'duration (seconds)')
print(paste("Coeficiente de variación: ", cv))

## [1] "Coeficiente de variación:  5.4511881438682"

latitude

varianza <- var(datos_us$latitude)
print(paste("Cuasivarianza: ", varianza))

## [1] "Cuasivarianza:  29.2851452840036"

sd <- sd(datos_us$latitude)
print(paste("Cuasidesviación típica: ", sd))

## [1] "Cuasidesviación típica:  5.41157512042506"

cv <- coef_var(datos_us$latitude)
print(paste("Coeficiente de variación: ", cv))

## [1] "Coeficiente de variación:  0.141580908541579"

longitude

varianza <- var(datos_us$longitude)
print(paste("Cuasivarianza: ", varianza))

## [1] "Cuasivarianza:  326.357120221342"

sd <- sd(datos_us$longitude)
print(paste("Cuasidesviación típica: ", sd))

## [1] "Cuasidesviación típica:  18.0653569082192"

cv <- coef_var(datos_us$longitude)
print(paste("Coeficiente de variación: ", cv))

## [1] "Coeficiente de variación:  -0.187626786221631"

datetime

Como no se pueden dividir fechas, no calculamos el Coeficiente de variación para datetime.

varianza <- var(datos_us$datetime)
print(paste("Cuasivarianza: ", varianza))

## [1] "Cuasivarianza:  108644846179268224"

sd <- sd(datos_us$datetime)
print(paste("Cuasidesviación típica: ", sd))

## [1] "Cuasidesviación típica:  329613176.586235"

Medidas de Forma

latitude

Al calcular la asimetría y curtosis de la latitud determinamos que: Obtenemos una asimetría negativa, por tanto la mayoría de datos se encuentran a la derecha de la media. Obtenemos una curtosis mayor que cero, por tanto, hay una mayor concentración de datos alrededor de la media

library(moments)
x <- (datos_us$latitude)
print(kurtosis(x))

## [1] 3.553226

print(skewness(x))

## [1] -0.04729688

print(mean(x))

## [1] 38.22249

hist(x)

### longitude Al calcular la asimetría y curtosis de la latitud determinamos que: Obtenemos una asimetría negativa, por tanto la mayoría de datos se encuentran a la derecha de la media. Obtenemos una curtosis mayor que cero, por tanto, hay una mayor concentración de datos alrededor de la media

library(moments)
x <- (datos_us$longitude)
print(kurtosis(x))

## [1] 2.315487

print(skewness(x))

## [1] -0.4947675

print(mean(x))

## [1] -96.28346

hist(x)

### datetime Al calcular la asimetría y curtosis de la latitud determinamos que: Obtenemos una asimetría negativa, por tanto la mayoría de datos se encuentran a la derecha de la media. Obtenemos una curtosis mayor que cero, por tanto, hay una mayor concentración de datos alrededor de la media

library(moments)
x <- (datos_us$datetime)
print(kurtosis(x))

## [1] 10.29702

print(skewness(x))

## [1] -2.502833

hist(x, 50)

### duration (seconds) Al calcular la asimetría y curtosis de la latitud determinamos que: Obtenemos una asimetría positiva, por tanto la mayoría de datos se encuentran a la izquierda de la media. Obtenemos una curtosis mayor que cero, por tanto, hay una mayor concentración de datos alrededor de la media

library(moments)
x <- (datos_us$'duration (seconds)')
print(kurtosis(x))

## [1] 264.1936

print(skewness(x))

## [1] 15.34129

print(mean(x))

## [1] 1086.33

Gráficos

Asociados a la Tabla de Frecuencias para algunos datos (excepto histogramas, usados anteriormente para ver mejor las medidas de forma) ### latitud

breaks <- seq(0, 90, by=10)
lat.cumfreq <- c(0,cumsum(table(cut(datos_us$latitude, breaks, right = FALSE))))

plot(breaks, lat.cumfreq,
     main="Latitud de los ovnis",
     xlab = "Latitud",
     ylab = "Latitudes acumuladas")
lines(breaks, lat.cumfreq)

### longitude

breaks <- seq(-200, 100, by=10)
longitude.cumfreq <- c(0,cumsum(table(cut(datos_us$longitude, breaks, right = FALSE))))
plot(breaks, longitude.cumfreq,
     main="Longitud de los ovnis",
     xlab = "Longitud",
     ylab = "Longitudes acumuladas")
lines(breaks, longitude.cumfreq)

Tema 4 - Probabilidad

Si quieres tratar de avistar un UFO la siguiente información que te proporcionamos te será de suma utilidad, como el pais, estado y ciudad donde debes ir para realizar un avistamiento, además de la forma que debes buscar y las horas a las que debes buscarlo, además de otros datos curiosos como la duración y los años que fueron más probables de avistar.

print("Probabilidad de avistar un ovni por país:")

## [1] "Probabilidad de avistar un ovni por país:"

dat <- read.csv("scrubbed.csv")
dat_country<-aggregate(dat$country, by=list(Country = dat$country), FUN=length)
dat$country <- ifelse(dat$country=="","otros",dat$country)
numAvistamientos<-sum(dat_country$x)
proobau<-dat_country[which(dat_country$Country == 'au'),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("Australia: ",probbauval,"%"))

## [1] "Australia:  0.669720659264054 %"

proobau<-dat_country[which(dat_country$Country == 'ca'),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("Canadá: ",probbauval,"%"))

## [1] "Canadá:  3.73450181745755 %"

proobau<-dat_country[which(dat_country$Country == 'de'),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("Alemania: ",probbauval,"%"))

## [1] "Alemania:  0.130707563611014 %"

proobau<-dat_country[which(dat_country$Country == 'gb'),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("Gran Bretaña: ",probbauval,"%"))

## [1] "Gran Bretaña:  2.37140865408555 %"

proobau<-dat_country[which(dat_country$Country == 'us'),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("Estados Unidos: ",probbauval,"%"))

## [1] "Estados Unidos:  81.056117113977 %"

proobau<-dat_country[which(dat_country$Country == ''),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("otros países: ",probbauval,"%"))

## [1] "otros países:  12.0375441916048 %"

maxcountry<-sum(by_country$x)

ggplot(by_country, aes(x=Country, y=(x/maxcountry)*100, fill=Country)) +
  geom_bar(stat="identity", width=1) +
  xlab('Paises') +
  ylab('Probabilidad') +
  guides(fill = guide_legend(title = "Country" ))

Debido a que es en Estados Unidos donde hay mayor probabilidad de avistar un UFO vamos a calcular los estados donde es mas probable lograrlo

datos_us$state <- ifelse(datos_us$state=="","otros",datos_us$state)

by_state <- aggregate(datos_us$state, by=list(State = datos_us$state), FUN=length)

avistTot<-sum(by_state$x)
mostseekstate<-by_state[which(by_state$x>16),]

maxstate<-by_state[which(by_state$x == max(by_state$x)),]
print(paste("El estado con mayor probabilidad de avistar un ovni es",maxstate$State,"con una probabilidad de un ",(maxstate$x/avistTot)*100,"%"))

## [1] "El estado con mayor probabilidad de avistar un ovni es ca con una probabilidad de un  13.8 %"

print(paste("A continuación se muestra una gráfica con los estados donde mayor probabilidad hay de avistar un ovni"))

## [1] "A continuación se muestra una gráfica con los estados donde mayor probabilidad hay de avistar un ovni"

ggplot(mostseekstate, aes(x=mostseekstate$State, y=(mostseekstate$x/avistTot)*100, fill=mostseekstate$State)) +
  geom_bar(stat="identity", width=1) +
  xlab('') +
  ylab('Probabilidad') +
  theme(axis.text.x=element_blank())+
  guides(fill = guide_legend(title = "Estados" ))

Ya que hemos calculado los estados haremos lo mismo con las ciudades

by_city <- aggregate(datos_us$city, by=list(City = datos_us$city), FUN=length)
avistTot<-sum(by_city$x)
mostseek<-by_city[which(by_city$x>4),]

maxcity<-by_city[which(by_city$x == max(by_city$x)),]
print(paste("La ciudad con mayor probabilidad de avistar un ovni es",maxcity$City,"con una probabilidad de un ",(maxcity$x/avistTot)*100,"%"))

## [1] "La ciudad con mayor probabilidad de avistar un ovni es phoenix con una probabilidad de un  1.2 %"

print(paste("A continuación se muestra una gráfica con las ciudades donde mayor probabilidad hay de avistar un ovni"))

## [1] "A continuación se muestra una gráfica con las ciudades donde mayor probabilidad hay de avistar un ovni"

ggplot(mostseek, aes(x=City, y=(x/avistTot)*100, fill=City)) +
  geom_bar(stat="identity", width=1) +
  xlab('') +
  ylab('Probabilidad') +
  theme(axis.text.x=element_blank())+
  guides(fill = guide_legend(title = "Ciudades" ))

En cuanto a la forma que tienes que buscar en el cielo para lograr un avistamiento exitoso

by_shape <- aggregate(datos_us$shape, by=list(shape = datos_us$shape), FUN=length)
by_shape$shape <- ifelse(by_shape$shape=="","otros",by_shape$shape)

max_shape<-sum(by_shape$x)
maxshape<-by_shape[which(by_shape$x == max(by_shape$x)),]
minshape<-by_shape[which(by_shape$x == min(by_shape$x)),]
print(paste("El ovni mas probabilidad de ser avistado son aquellos ovnis con forma",maxshape$shape,"con una probabilidad del",(maxshape$x/avistTot)*100,"%"))

## [1] "El ovni mas probabilidad de ser avistado son aquellos ovnis con forma light con una probabilidad del 21.2 %"

print(paste("Mientras que los más raros con mayor dificultad de avistamiento son aquellos con forma"))

## [1] "Mientras que los más raros con mayor dificultad de avistamiento son aquellos con forma"

print (paste(minshape$shape))

## [1] "cross"

print(paste("con una probabilidad de avistamiento del",(minshape$x[1]/avistTot)*100,"%"))

## [1] "con una probabilidad de avistamiento del 0.2 %"

ggplot(by_shape, aes(x=by_shape$shape,y=(by_shape$x/max_shape)*100,fill=by_shape$shape)) +
  geom_bar(stat="identity", width=1) +
  xlab('') +
  ylab('Avistamientos') +
theme(axis.text.x=element_blank())+
  guides(fill = guide_legend(title = "Formas" ))

La hora también es importante, hay que saber cuando buscar y cuando descansar

by_horas<- aggregate(datos_us$datetime, by=list(Hour = format(datos_us$datetime,"%H")), FUN=length)
maxhora<-by_horas[which(by_horas$x == max(by_horas$x)),]
maxhour<-sum(by_horas$x)
print(paste("La hora en la que es más probable avistar un ovni son las",maxhora$Hour,':00',"con una probabilidad del",(maxhora$x/avistTot)*100,"%"))

## [1] "La hora en la que es más probable avistar un ovni son las 22 :00 con una probabilidad del 14.4 %"

plot(by_horas$Hour,(by_horas$x/maxhour)*100,type="s",col="dark red", xlab = "Horas del día", ylab = "Probabilidad", main = "Probabilidad de avistamiento por hora")

También es interesante ver cuando hubo una mayor probabilidad de avistar un ovni

anios<- aggregate(datos_us$datetime, by=list(Date = format(datos_us$datetime,"%Y")), FUN=length)

maxanio<-anios[which(anios$x == max(anios$x)),]
maxyear<-sum(anios$x)
print(paste("El año con mayor probabilidad en caso de avistamiento fué",maxanio$Date,"con una probabilidad del",(maxanio$x/avistTot)*100,"%"))

## [1] "El año con mayor probabilidad en caso de avistamiento fué 2012 con una probabilidad del 9.9 %"

plot(anios$Date,(anios$x/maxyear)*100,type="l",col="dark blue", xlab = "Linea temporal", ylab = "Probabilidad", main = "Probabilidad de avistamientos por año")

Observamos la probabilidad de lograr un avistamiento de larga duración

durability<- aggregate(datos_us$`duration (seconds)`, by=list(duracion = datos_us$`duration (seconds)`), FUN=length)
  maxduracion<-sum(durability$x)  
  maxdur<-durability[which(durability$x>1),]
  
durationProb<-durability[which(durability$x<=10),]
sumaprobduracion<-sum(durationProb$x)
print(paste("La probabilidad de ver un UFO durante menos de 10 segundos es del",(sumaprobduracion/maxduracion)*100,"%"))

## [1] "La probabilidad de ver un UFO durante menos de 10 segundos es del 12.3 %"

durationProb<-durability[which(durability$x<=60),]
sumaprobduracion<-sum(durationProb$x)
print(paste("La probabilidad de ver un UFO durante menos de un minuto es del",(sumaprobduracion/maxduracion)*100,"%"))

## [1] "La probabilidad de ver un UFO durante menos de un minuto es del 60.9 %"

durationProb<-durability[which(durability$x>60),]
sumaprobduracion<-sum(durationProb$x)
print(paste("La probabilidad de ver un UFO durante mas de un minuto es del",(sumaprobduracion/maxduracion)*100,"%"))

## [1] "La probabilidad de ver un UFO durante mas de un minuto es del 39.1 %"

plot(maxdur$duracion,(maxdur$x/sumaprobduracion)*100,type="l",col="blue", xlab = "Segundos", ylab = "Probabilidad", main = "Probabilidad de avistamientos por más de un segundo")

A continuación realizaremos algunas operaciones con sucesos de probabilidad

prob1 <- datos_us[which(datos_us$city == 'seattle'&format(datos_us$datetime,"%H")>12),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en seattle por la tarde",(sum(prob2$x)/avistTot)*100,"%"))

## [1] "La probabilidad de ver un UFO en seattle por la tarde 0.5 %"

prob1 <- datos_us[which(datos_us$city == 'seattle'|format(datos_us$datetime,"%H")==21),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en seattle o verlo a las 21:00 es del",(sum(prob2$x)/avistTot)*100,"%"))

## [1] "La probabilidad de ver un UFO en seattle o verlo a las 21:00 es del 13.7 %"

prob1 <- datos_us[which(datos_us$city != 'seattle'&format(datos_us$datetime,"%H")!=21),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en una ciudad que no sea seattle a una hora distinta a las 21:00 es del",(sum(prob2$x)/avistTot)*100,"%"))

## [1] "La probabilidad de ver un UFO en una ciudad que no sea seattle a una hora distinta a las 21:00 es del 86.3 %"

prob1 <- datos_us[which(datos_us$city != 'seattle'&format(datos_us$datetime,"%H")==21),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en una ciudad que no sea seattle a las 21:00 es del",(sum(prob2$x)/avistTot)*100,"%"))

## [1] "La probabilidad de ver un UFO en una ciudad que no sea seattle a las 21:00 es del 13.2 %"

prob1 <- datos_us[which(datos_us$city != 'seattle'&format(datos_us$datetime,"%H")==21&(datos_us$shape == 'light'|datos_us$shape=='unknown')),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en una ciudad distinta de seattle a las 21:00 con una forma que sea lumínica o desconocida",(sum(prob2$x)/avistTot)*100,"%"))

## [1] "La probabilidad de ver un UFO en una ciudad distinta de seattle a las 21:00 con una forma que sea lumínica o desconocida 4 %"

prob1 <- datos_us[which(datos_us$city == 'seattle'&format(datos_us$datetime,"%H")!=21&(datos_us$shape != 'light'&datos_us$shape!='unknown')&datos_us$`duration (seconds)`>10),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en seattle a una hora distinta de las 21:00 con una forma que sea distinta de lumínica o desconocida pero con una duracion de mas de 10 segundos",(sum(prob2$x)/avistTot)*100,"%"))

## [1] "La probabilidad de ver un UFO en seattle a una hora distinta de las 21:00 con una forma que sea distinta de lumínica o desconocida pero con una duracion de mas de 10 segundos 0.3 %"

Respondemos algunas preguntas de probabilidad condicional

maxdat<-count(datos)
prob1 <- datos_us[which(format(datos_us$datetime,"%H")==15),]
prob2<-count(prob1)
print(paste("Dado que son las 15:00 que probabilidad hay de ver un UFO?",(prob2/avistTot)*100,"%"))

## [1] "Dado que son las 15:00 que probabilidad hay de ver un UFO? 1.6 %"

prob1 <- datos_us[which(datos_us$city == 'richmond'),]
prob2<-count(prob1)
print(paste("Si vivo en Richmond que probabilidad tengo de ver un UFO?",(prob2/avistTot)*100,"%"))

## [1] "Si vivo en Richmond que probabilidad tengo de ver un UFO? 0.3 %"

prob1 <- datos_us[which(datos_us$shape=='changing'),]
prob2<-count(prob1)
print(paste("En caso de ver un ovni que probabilidad hay de que sea cambiante?",(prob2/avistTot)*100,"%"))

## [1] "En caso de ver un ovni que probabilidad hay de que sea cambiante? 2.5 %"

prob1 <- datos[which(datos$country != 'us'),]
prob2<-count(prob1)
print(paste("Si vivo fuera de Estados Unidos que probabilidad hay de avistar un ovni?",(prob2/maxdat)*100,"%"))

## [1] "Si vivo fuera de Estados Unidos que probabilidad hay de avistar un ovni? 18.9425722359854 %"

Tema 5 - Variables Aleatorias y Modelos de Probabilidad

Distribucion normal

Como los datos, excepto los extremos, están aproximadamente sobre la linea, podemos decir que exluyendo los extremos, los datos siguen una distribucion normal.

dia <-as.numeric(format(datos_us$datetime,'%d'))

qqnorm(dia, pch = 1, frame = FALSE)
qqline(dia, col = "steelblue", lwd = 2)

También tenemos la media, moda, varianza y desviación típica de estos datos.

mean(dia)

## [1] 14.984

mode(dia)

## [1] 15

var(dia)

## [1] 79.80155

sd(dia)

## [1] 8.933171

Funcion de distribucion acumulada respecto a los días de los mes mas comunes de todos los años en el que se han visto ovnis en EEUU.

f1 <- pnorm(dia,mean(dia),sd(dia))

plot(dia,f1,main = "Funcion de distribucion acumulada")

Calculamos la probabilidad de que se avisten ovnis durante los primeros 10 dias de cada mes

P(X≤10)

pnorm(10,mean(dia),sd(dia))

## [1] 0.2884493

Calculamos la probabilidad de que se avisten ovnis durante los ultimos 5 dias de cada mes

P(X>25)

1-pnorm(25,mean(dia),sd(dia))

## [1] 0.1310983

Calculamos la probabilidad de que se avisten ovnis entre los dias 10 y 20 de cada mes.

P(10≤X≤20)

pnorm(20,mean(dia),sd(dia)) - pnorm(10,mean(dia),sd(dia))

## [1] 0.4243233

Podemos representar este último gráficamente:

regionX=seq(10,20,0.01)
xP <- c(10,regionX,20)
yP <- c(0,dnorm(regionX,15,12),0)

curve(dnorm(x,15,12),xlim=c(0,30),yaxs="i",ylim=c(0,0.035),ylab="f(x)",
      main='Densidad P(10<X<20)') 
polygon(xP,yP,col="orange1")
box()

Diagrama de Poisson

Podemos averiguar cuantos ovnis fueron divisados en 2012.

year <-as.numeric(format(datos_us$datetime,'%y'))

landa <- sum(year == 12)

plot(dpois(0:200, landa), type = "h", lwd = 2,
     main = "Función de masa de probabilidad",
     ylab = "P(X = x)", xlab = "Número de ovnis avistados en 2012")

Distribución Exponencial

Una vez calculada la tasa de llegada de los ovnis en 2012 con la distribucion de Poisson, podemos dibujar la gráfica de la función de densidad exponencial.

# Rejilla del eje X
x <- seq(0, 0.1, 0.01)

dexp(x,landa)

##  [1] 99.000000000 36.786092411 13.668854494  5.079027723  1.887248315
##  [6]  0.701257484  0.260570935  0.096822086  0.035976830  0.013368152
## [11]  0.004967294

# lambda
plot(x, dexp(x, landa), type = "l",
     ylab = "f(x)", lwd = 2, col = "red")

De la misma forma tenemos la funcion de distribución de exponencial acumulada

# Rejilla de valores del eje X
x <- seq(0,0.1, 0.01)

pexp(x,landa)

##  [1] 0.0000000 0.6284233 0.8619308 0.9486967 0.9809369 0.9929166 0.9973680
##  [8] 0.9990220 0.9996366 0.9998650 0.9999498

# lambda
plot(x, pexp(x, landa), type = "l",
     ylab = "F(x)", lwd = 2, col = "orange")

La probabilidad de que se divisen menos de 100 ovnis en 2012 P(X<=100) es:

ppois(100,landa)

## [1] 0.5663565

Probabilidad de que se divisen entre 80 y 95 ovnis en 2012 P(80< X< 95)

A modo ilustrativo, podemos representarlo en una gráfica

ppois(95,landa) - ppois(80,landa)

## [1] 0.3396651

pois_sum <- function(lambda, lb, ub, col = 4, lwd = 1, ...) {
    x <- 0:(lambda + lambda * 2)
    
    if (missing(lb)) {
       lb <- min(x)
    }
    if (missing(ub)) {
        ub <- max(x)
    }
      
    plot(dpois(x, lambda = lambda), type = "h", lwd = lwd, ...)
  
    if(lb == min(x) & ub == max(x)) {
        color <- col
    } else {
        color <- rep(1, length(x))
        color[(lb + 1):ub ] <- col
    }
    
    lines(dpois(x, lambda = lambda), type = "h",
          col =  color, lwd = lwd, ...)
}

pois_sum(landa,80, 95, lwd = 2,
           col = 2, ylab = "P(X = x)", xlab = "Ovnis avistados en 2012")

Gráfico de la función de distribución exponencial

plot(ppois(40:100, landa), type = "s", lwd = 2,
     main = "Función de distribución exponencial",
     xlab = "Número de eventos", ylab = "F(x)")

También podemos obtener los cuantiles correspondientes de la distribucion exponencial

plot(qpois(seq(0,1,0.001), landa),
     main = "Función cuantil",
     ylab = "Q(p)", xlab = "Cuantiles",
     type = "s", col = 3, xaxt = "n")

axis(1, labels = seq(0, 1, 0.1), at = 0:10 * 100)

Distribucion Binomial

Para los datos, vamos a calcular la probabilidad de se vean ovnis el dia 24 de cada mes.

Con lo que si se ven ovnis los días 24, son aciertos, y al contrario son fallos.

dia <-as.numeric(format(datos_us$datetime,'%d'))

tam<- 1000
x <- (dia == 24)

y <- sum(x)
  
probEnsayo <- (y/tam)

print(probEnsayo)

## [1] 0.03

Funcion de probabilidad binomial

plot(dbinom(1:100, tam, probEnsayo), type = "h", lwd = 2,
     main = "Función de probabilidad binomial",
     ylab = "P(X = x)", xlab = "Número de éxitos")

Probabilidad de encontrar a lo largo de 69 años (2014 - 1945) mas de 35 ovnis el día 24 de cada mes:

P(X > 35)

1 - pbinom(35, tam, probEnsayo)

## [1] 0.1539237

Probabilidad de encontrar a lo largo de 69 años (2014 - 1945) menos de 25 ovnis el día 24 de cada mes:

P(X <= 25)

pbinom(25, tam, probEnsayo)

## [1] 0.2044652

Ahora graficmente la probabilidadd de encontrar entre 25 y 35 ovnis.

P(25 <= X <= 35)

binom_sum <- function(size, prob, lb, ub, col = 4, lwd = 1, ...) {
    x <- 0:size
    
    if (missing(lb)) {
       lb <- min(x)
    }
    if (missing(ub)) {
        ub <- max(x)
    }
      
    plot(dbinom(x, size = size, prob = prob), type = "h", lwd = lwd, ...)
  
    if(lb == min(x) & ub == max(x)) {
        color <- col
    } else {
        color <- rep(1, length(x))
        color[(lb + 1):ub ] <- col
    }
    
    lines(dbinom(x, size = size, prob = prob), type = "h",
          col =  color, lwd = lwd, ...)
}

binom_sum(0 :100, probEnsayo, lb = 25, ub = 35, lwd = 2,
          ylab = "P(X = x)", xlab = "Número de éxitos")

Funcion de distribucion binomial

plot(pbinom(10:50, tam, probEnsayo), type = "s", lwd = 2,
     main = "Función de distribución binomial",
     xlab = "Número de éxitos", ylab = "F(x)")

Teorema central del límite

Hemos usado la distribucion binomial para el teorema central del límite,

Realizamos el histograma de las 6 últimas muestras y de la media muestral,

tamMuestra=1000
numMuestras=25
#hemo 

#Se deja ese hueco vacío aunque salga warning
mat=matrix(, nrow = numMuestras, ncol = tamMuestra)
mediax=vector()

# Genera numMuestras muestras aleatorias 
# Para cada muestra, almacena la media en el vector anterior
for (i in 1:numMuestras){
  x=rbinom(tamMuestra,200,probEnsayo)
  mat[i,]=x
  mediax[i]=mean(x)
}

# Muestra la media de las 6 Ãºltimas muestras generadas
print(c(mediax[numMuestras],mediax[numMuestras-1],mediax[numMuestras-2],
        mediax[numMuestras-3],mediax[numMuestras-4],mediax[numMuestras-5]))

## [1] 6.038 6.027 5.995 5.961 5.875 5.935

# Muestra el histograma de las 6 Ãºltimas muestras y de la media muestral
par(mfrow=c(3, 3))
hist(mat[numMuestras,],probability=TRUE)
hist(mat[numMuestras-1,],probability=TRUE)
hist(mat[numMuestras-2,],probability=TRUE)
hist(mat[numMuestras-3,],probability=TRUE)
hist(mat[numMuestras-4,],probability=TRUE)
hist(mat[numMuestras-5,],probability=TRUE)
hist(mediax,probability=TRUE)

Tema 6 - Muestreo, Estimación Puntual y por Intervalos de Confianza

Vamos a construir un intervalo de confianza del 95% para la duracion media de los avistamientos en EEUU.

#filtramos por duracion menor o igual a 7200, para descartar valores muy altos que son inútiles en nuestro análisis
retval <- subset(datos_us, datos_us$`duration (seconds)` < 7201, select = c(state, `duration (seconds)`))

intervalo <- t.test(retval$'duration (seconds)', conf.level = 0.95)#sacamos el intervalo de confianza
print(intervalo$conf.int)#intervalo de confianza

## [1] 490.9120 624.0471
## attr(,"conf.level")
## [1] 0.95

print(mean(retval$'duration (seconds)'))#media

## [1] 557.4796

print(paste("El resultado indica que la media de la variable duracion (seconds) es de ", mean(retval$duration..seconds.), " el cual se encuentra con una confianza del 95% en el intervalo ", intervalo$conf.int[1], " ", intervalo$conf.int[2]))

## [1] "El resultado indica que la media de la variable duracion (seconds) es de  NA  el cual se encuentra con una confianza del 95% en el intervalo  490.912046871598   624.047136801871"

car::qqPlot(retval$'duration (seconds)', pch=19,
       main='QQplot para la duración de avistamientos',
       xlab='Cuantiles teóricos',
       ylab='Cuantiles muestrales')

## [1]  99 140

hist(retval$'duration (seconds)', freq=TRUE,
     main='Histograma para la duración de avistamientos',
     xlab='Duracion (s)',
     ylab='Frecuencia')

     #Podemos observar que la mayoría de avistamientos duran muy poco tiempo.

Tema 7 - Contrastes de Hipótesis

Dado que la mayoría de los datos son de EEUU, queremos saber si los OVNIs realmente prefieren estar ahí. Para esto, haremos un contraste de la duración del avistamiento entre EEUU y todos los demás. Podemos asumir que como los datos son del National UFO Reporting Center (NUFORC) basado EEUU, datos de otros países no son datos poblacionales.

Basaremos el estudio en las siguentes hipótesis:

Hipótesis	Contraste
H0 (Nula)	µ(EEUU) = µ(Otros)
H1 (Alternativa)	µ(EEUU) != µ(Otros)

Donde µ(EEUU) y µ(Otros) son la duración media de un avistamiento en EEUU y en otros países respectivamente. Tomaremos el nivel de significancia α = 0.1

alpha <- 0.1

datos_us <- datos[which(datos$country == 'us'),]
datos_us <- datos_us %>% select(-country)

datos_otros <- datos[which(datos$country != 'us'),]
datos_otros <- datos_otros %>% select(-country)

# Eliminamos valores extremos
duracion_us <- datos_us$`duration (seconds)`[!(datos_us$`duration (seconds)` %in% boxplot.stats(datos_us$`duration (seconds)`)$out)]
duracion_otros <- datos_otros$`duration (seconds)`[!(datos_otros$`duration (seconds)` %in% boxplot.stats(datos_otros$`duration (seconds)`)$out)]

Como tenemos un número de datos muy grande, usamos el teorema central del límite (verificando que es distribución normal) para calcular µ(EEUU). Dado que las duraciones son muy variables, asumimos que las desviaciones desconocidas son distintas.

muestra_duracion_us <- replicate(100, sample(duracion_us, size = 1000, replace = TRUE) %>% mean())

car::qqPlot(muestra_duracion_us)

## [1] 70 97

var_us <- var(muestra_duracion_us)
x_us <- mean(muestra_duracion_us)

Hacemos lo mismo para las duraciones de otros países

muestra_duracion_otros <- replicate(100, sample(duracion_otros, size = 1000, replace = TRUE) %>% mean())

car::qqPlot(muestra_duracion_otros)

## [1] 96 92

var_otros <- var(muestra_duracion_otros)
x_otros <- mean(muestra_duracion_otros)

Calculamos la region de rechazo

f <- ((var_us / 100 + var_otros / 100) ^ 2) / ((var_us/100)^2/99 + (var_otros/100)^2/99)
f <- round(f)
cuantil <- qt(c(.025, .975), df=f)

print(paste("La region de rechazo es T < ", cuantil[1], " o T > ", cuantil[2]))

## [1] "La region de rechazo es T <  -1.97201747783631  o T >  1.97201747783631"

estadistico <- abs(x_us - x_otros) / sqrt(var_us / 100 + var_otros / 100)
print(estadistico)

## [1] 5.945632

Como nuestro estadistico está en la region de rechazo, rechazamos la hipótesis nula. Eso significa que efectivamente, la duración de un avistamiento es mayor en EEUU.

Podemos verificar esto usando el test t de student ya que nuestra n > 30

t.test(x=muestra_duracion_us, y=muestra_duracion_otros, alternative = "two.sided", conf.level  = 0.9)

## 
##  Welch Two Sample t-test
## 
## data:  muestra_duracion_us and muestra_duracion_otros
## t = 5.9456, df = 197.71, p-value = 1.228e-08
## alternative hypothesis: true difference in means is not equal to 0
## 90 percent confidence interval:
##   6.158525 10.899955
## sample estimates:
## mean of x mean of y 
##  264.8760  256.3468

EA - Analisis OVNIs

Objetivos

Datos

Tema 3 - Estadística Descriptiva y Regresión

Medidas de Posición - Tendencia Central y No Central

datetime

latitude

longitude

duration (seconds)

Medidas de Dispersión

duration (seconds)

latitude

longitude

datetime

Medidas de Forma

latitude

Gráficos

Tema 4 - Probabilidad

Tema 5 - Variables Aleatorias y Modelos de Probabilidad

Distribucion normal

Diagrama de Poisson

Distribución Exponencial

Distribucion Binomial

Teorema central del límite

Tema 6 - Muestreo, Estimación Puntual y por Intervalos de Confianza

Tema 7 - Contrastes de Hipótesis